rm(list=ls()) # clear workspace

# set working directory
setwd("C:/Users/cailin.slattery/Dropbox/Incentives/postJM/Replication")

set.seed(31389) 
library(mvtnorm)
library(MASS)

######################### DEFINE FUNCTIONS #############################
maxn <- function(n) function(x) order(x, decreasing = TRUE)[n]

process_sim <- function(pi, v, N2) {
   w = v+pi
   #max welfare
   max_w = apply(w,1, max)
   w_loc = cbind(matrix(1:N2, N2,1),apply(w,1,which.max))
   # profit in max welfare place
   win_pi= pi[w_loc]
   # value in max welfare place
   win_v= v[w_loc]
   
   #welfare in max pi place
   pi_loc = cbind(matrix(1:N2, N2,1),apply(pi,1,which.max))
   max_piw= w[pi_loc]
   max_pi= pi[pi_loc]
   v0 = v[pi_loc]
   
   #change in welfare
   delta_w = max_w - max_piw
   
   #now figure out split between firm and state
   #second highest welfare (runner-up)
   ru_loc = cbind(matrix(1:N2, N2,1),apply(w,1, maxn(2)))
   # value in runner-up place
   ru_v = v[ru_loc]
   # profit in runner-up place
   ru_pi = pi[ru_loc]
   
   subsidy = ru_pi - win_pi + ru_v 
   firm_payoff = subsidy + win_pi
   state_payoff = win_v - subsidy
   delta_firm = firm_payoff - max_pi
   delta_state = state_payoff - v0
   
   results_base = data.frame(subsidy, delta_w, max_w, max_piw, max_pi, win_pi, v0, win_v)
}

############################################################################################################### 
# PRIMITIVES
############################################################################################################### 
v_data <- read.csv("data/Hv_manuf.csv")
fit_v <- fitdistr(v_data$x, "lognormal")
pi_data <- read.csv("data/manufpi.csv")
fit_pi <- fitdistr(pi_data$x, "normal")

###############################################################################################################
###############################################################################################################
#################BASELINE: Change variance of v
###############################################################################################################
###############################################################################################################

N=10
r <-seq(from = 0 , to = 2*fit_v$estimate[2], by = .1)
N2=length(r)*1000
testv = lapply(r, function(x) {matrix(rlnorm(N*1000, fit_v$estimate[1], x),ncol=N)})
v <-do.call(rbind,testv)
v[which(v>1500)]<-mean(v)
pi1 = matrix(rnorm(N*1000,fit_pi$estimate[1],fit_pi$estimate[2]), ncol=N)
pi1[which(pi1<0)]<-0
pi <- do.call(rbind, replicate(N2/1000, pi1, simplify=FALSE))

results_base <- process_sim(pi,v, N2)

#export results from welfare simulation
write.csv(results_base, "data/sim/results_varv10.csv")

N=5
testv = lapply(r, function(x) {matrix(rlnorm(N*1000, fit_v$estimate[1], x),ncol=N)})
v <-do.call(rbind,testv)
v[which(v>1500)]<-mean(v)
pi1 = matrix(rnorm(N*1000,fit_pi$estimate[1],fit_pi$estimate[2]), ncol=N)
pi1[which(pi1<0)]<-0
pi <- do.call(rbind, replicate(N2/1000, pi1, simplify=FALSE))

results_base <- process_sim(pi,v, N2)

#export results from welfare simulation
write.csv(results_base, "data/sim/results_varv5.csv")

N=3
testv = lapply(r, function(x) {matrix(rlnorm(N*1000, fit_v$estimate[1], x),ncol=N)})
v <-do.call(rbind,testv)
v[which(v>1500)]<-mean(v)
pi1 = matrix(rnorm(N*1000,fit_pi$estimate[1],fit_pi$estimate[2]), ncol=N)
pi1[which(pi1<0)]<-0
pi <- do.call(rbind, replicate(N2/1000, pi1, simplify=FALSE))

results_base <- process_sim(pi,v, N2)

#export results from welfare simulation
write.csv(results_base, "data/sim/results_varv3.csv")

N=2
testv = lapply(r, function(x) {matrix(rlnorm(N*1000, fit_v$estimate[1], x),ncol=N)})
v <-do.call(rbind,testv)
v[which(v>1500)]<-mean(v)
pi1 = matrix(rnorm(N*1000,fit_pi$estimate[1],fit_pi$estimate[2]), ncol=N)
pi1[which(pi1<0)]<-0
pi <- do.call(rbind, replicate(N2/1000, pi1, simplify=FALSE))

results_base <- process_sim(pi,v, N2)

#export results from welfare simulation
write.csv(results_base, "data/sim/results_varv2.csv")

###############################################################################################################
###############################################################################################################
###################### Change Correlation between v and pi
###############################################################################################################
###############################################################################################################

N=10
r <- seq(from =-0.5, to = 0.5, by = .02)
N2=length(r)*1000

s_v <- log(fit_v$estimate[2])
s_pi <- fit_pi$estimate[2]

s <- lapply(r, function(x) {matrix(c(s_v, x*sqrt(s_v)*sqrt(s_pi),  x*sqrt(s_v)*sqrt(s_pi), s_pi), nrow = 2)})
x <- lapply(s, function(x) {rmvnorm(N*1000, mean = c(fit_v$estimate[1],fit_pi$estimate[1]), sigma = x, method = "chol")})
xall<- do.call(rbind, x)
v <- matrix(exp(xall[,1]), nrow=N2, byrow=TRUE)
pi <- matrix(xall[,2], nrow=N2, byrow=TRUE)
v[which(v>1500)]<-mean(v)
pi[which(pi<0)]<-0

results_base <- process_sim(pi,v, N2)

#export results from welfare simulation
write.csv(results_base, "data/sim/results_corr10.csv")

N=5
s <- lapply(r, function(x) {matrix(c(s_v, x*sqrt(s_v)*sqrt(s_pi),  x*sqrt(s_v)*sqrt(s_pi), s_pi), nrow = 2)})
x <- lapply(s, function(x) {rmvnorm(N*1000, mean = c(fit_v$estimate[1],fit_pi$estimate[1]), sigma = x, method = "chol")})
xall<- do.call(rbind, x)
v <- matrix(exp(xall[,1]), nrow=N2, byrow=TRUE)
pi <- matrix(xall[,2], nrow=N2, byrow=TRUE)
v[which(v>1500)]<-mean(v)
pi[which(pi<0)]<-0

results_base <- process_sim(pi,v, N2)

#export results from welfare simulation
write.csv(results_base, "data/sim/results_corr5.csv")

N=3
s <- lapply(r, function(x) {matrix(c(s_v, x*sqrt(s_v)*sqrt(s_pi),  x*sqrt(s_v)*sqrt(s_pi), s_pi), nrow = 2)})
x <- lapply(s, function(x) {rmvnorm(N*1000, mean = c(fit_v$estimate[1],fit_pi$estimate[1]), sigma = x, method = "chol")})
xall<- do.call(rbind, x)
v <- matrix(exp(xall[,1]), nrow=N2, byrow=TRUE)
pi <- matrix(xall[,2], nrow=N2, byrow=TRUE)
v[which(v>1500)]<-mean(v)
pi[which(pi<0)]<-0

results_base <- process_sim(pi,v, N2)

#export results from welfare simulation
write.csv(results_base, "data/sim/results_corr3.csv")

N=2
s <- lapply(r, function(x) {matrix(c(s_v, x*sqrt(s_v)*sqrt(s_pi),  x*sqrt(s_v)*sqrt(s_pi), s_pi), nrow = 2)})
x <- lapply(s, function(x) {rmvnorm(N*1000, mean = c(fit_v$estimate[1],fit_pi$estimate[1]), sigma = x, method = "chol")})
xall<- do.call(rbind, x)
v <- matrix(exp(xall[,1]), nrow=N2, byrow=TRUE)
pi <- matrix(xall[,2], nrow=N2, byrow=TRUE)
v[which(v>1500)]<-mean(v)
pi[which(pi<0)]<-0

results_base <- process_sim(pi,v, N2)

#export results from welfare simulation
write.csv(results_base, "data/sim/results_corr2.csv")



